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ABSTRACT 

The existence of a bullet cluster (such as 1E0657-56) poses a challenge to the 
concordance A cold dark matter model. Here we investigate the velocity distribu- 
tion of dark matter halo pairs in large A'^-body simulations with differing box sizes 
{250h^^ Mpc — 2ft,^^Gpc) and resolutions. We examine various basic statistics such as 
the halo masses, pairwise halo velocities (W12), coUisional angles, and pair separation 
distances. We then compare our results to the initial conditions required to repro- 
duce the observational properties of 1E0657-56 in non-cosmological hydrodynamical 
simulations. 

We find that the high velocity tail of the W12 distribution extends to greater veloc- 
ities as we increase the simulation box size. We also find that the number of high-wi2 
pairs increases as we increase the particle count and resolution with a fixed box size, 
however, this increase is mostly due to lower mass halos which do not match the 
observed masses of 1E0657-56. We find that the redshift evolution effect is not very 
strong for the i;i2 distribution function between 2=0.0 and ^'^0.5. 

We identify some pairs whose W12 resemble the required initial conditions, however, 
even the best candidates have either wrong halo mass ratios, or too large separations. 
Our simulations suggest that it is very difficult to produce such initial conditions at 
z = 0.0,0.296, & 0.489 in comoving volumes as large as {2h~^Gpc)^. Based on the 
extrapolation of our cumulative f 12 function, we find that one needs a simulation with 
a comoving box size of (4.48 ft-~^ Gpc)'^ and 2240"^ DM particles in order to produce 
at least one pair of halos that resembles the required W12 and observed masses of 
1E0657-56. From our simulated W12 probability distribution function, we find that 
the probability of finding a halo pair with i;i2 ^ 3000 km s~^ and masses ^ lO^^M© 
to be 2.76 x 10"^ at ^=0.489. We conclude that either 1E0657-56 is incompatible 
with the concordance ACDM universe, or the initial conditions suggested by the non- 
cosmological simulations must be revised to give a lower value of W12 . 

Key words: method : N-body simulations — galaxies : evolution — galaxies : for- 
mation — galaxies: clusters — cosmology : theory — cosmology : dark matter 



1 INTRODUCTION 

It is widely believed that the structure formation in our Uni- 
verse is largely driven by the gravity of dark matter. There- 
fore it is worthwhile to probe dark matter dynamics through 
measurements of galaxy peculiar velocities and constrain our 
cosmological model by comparing against numerical simu- 
lations. In fact there has been extensive work along these 
lines, recovering the local density field from the measured 
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velocity field (iBertschinger fc Dekelll 19891: iDavis et al.]|l996l : 



IWillick et al.lll996l ). Unfortunately the observations of pecu- 
liar velocity fields contain large uncertainties, and accurate 
determination of the cosmological mass density parameter 
Sim turned out to be difficult using this method. 

More recently, clusters of galaxies have been used to 
prove the existence of dark matter itself, thanks to accu- 
rate measurements of projected dark matter density using 
weak and strong lensing techniques. Some clusters show 
signs of a cluster-cluster merger, where the baryonic com- 
ponent and coUisionless dark matter show different spatial 
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distributions, strongly supporting the existence of dark mat- 
ter. Furthermore, using the shock features seen in the gas, 
one can infer the c ollision velocity of two ga laxy clusters 
l|Clowe et alj |2004 l2006l : iBradac et all |2006| ). These new 
observations have brought renewed interest to dark mat- 
ter dynamics and using it to check th e standard A cold dark 
matter (ACDM) cosmologic al model (jEfstathiou et al.ll 19901 : 
lOstriker fc Steinhardtlll995D . 

In particular, the observations of the massive cluster 
of galaxies 1E0657-56 seem to suggest a much higher rel- 
ative dark matter halo velocity than one would expect in 
the ACDM model. This system includes a massive sub- 
cluster (the "bullet") with Mbuiiot ^ 1.5 x lO^Af© that 
has fallen through the parent cluster of Mparcnt — 1.5 x 
W^^Mq roughly 150 million years ago, and is separated by 
~ 0.72 Mpc on the sky at an observed redsh ift of 2=0.296 
ICloweet al.ll2004 l2006l : iBradac et all [200^ ). The unique- 
ness of this system comes from the collision trajectory being 
almost perpendicular to our line of sight. This provides an 
opportunity to better study the dynamics of large cluster 
collisions. The Chandra observations revealed that the pri- 
mary baryonic component had been stripped away in the 
collision and resided betw een the two clust ers in the form 
of hot X-ray emitting gas (|Markevitchll2006l ) . This provides 
strong evidence for the existence of dark matter (DM); As 
the two clusters passed through each other, the baryonic 
components interacted and slowed down due to ram pres- 
sure, while the dark matter component was allowed to move 
ahead of the gas since it only interacts through gravity with- 
out dissipation. One can infer the velocity of the bow shock 
preceding the 'bullet' using the shock Mach number and 
a measurement of the pre-shock temperature. The inferred 
shock vel ocity was found to be iighock ~ 4740I55Q kms~^ 
(JMarkevitch 2006,). 



Havashi fc White! (|2006l ) examined the Millennium Run 

ISpringel et al.ll2005l ) in search for such a sub-cluste 



ter mov- 
ing with a velocity relative to its parent c luster of Wbuiict ~ 
4:500tllo° km s"^ (JMarkevitch et al.ll2004l ). Due to the lim- 
ited volume of the simulation (500 h~^Mpc)^, few halos had 
masses comparable to 1E0657-56. Still they estimated that 
about 1 in 100 have velocities comparable to the bullet clus- 
ter, and concluded that the event is not impossible within 
the current ACDM model. 

It is often assumed that the inferred shock velocity is 
equal to the velocity of the dark matter 'bullet' itself. Several 
groups have shown, however, that this is not necessarily true 
throu gh the use of non-cosmolo gical hydrodynamic simula- 
tions. |Mi]osaAdjevic_et_alJ (120071 ) used two dimensional simu- 
lations to find that the subcluster's velocity differed from the 
shock velocity by about 16%, bringing the relative velocity of 
DM halos down to ~ 3980 kms~^. They assumed a zero rela- 
tive velocity at a separation distance of 4.6 Mpc for their ini- 
tial conditions. They also emphasized that their conclusion 
is sensitive t o the initial mass and ga s density profile of the 
two clusters. ISpringel fc FarraJ (|2007l ) was able to reproduce 
the inferred shock velocity through the use of an idealized 
three dimensional hydrodynamic simulation with initial con- 
ditions that assumed a relative velocity of 2057 km s"'^ at a 
separation distance of 3.37 Mpc, and found that the subclus- 
ter was moving with a relative speed of only ^2 6 00 km s~^ 
just after the collision. Mastropietro fc Burkerd (|200a ) ar- 
gued that ISpringel fc Farran (|2007l ) failed to reproduce the 



observed displacement of X-ray peaks that represent an im- 
portant indicator of the nature of the interaction. In their 
simulations they found that in order to reproduce the ob- 
servational data of 1E0657-56 a relative halo infall velocity 
of ~3000kms~^ at an initial separation distance of 5 Mpc 
was required. 



E 



Similar to previou s work by iHavashi fc Whitd (|2006l ) , 
ee fc Komatsul (|20ld ) quantified the likelihood of finding 



bullet-like s ystems in the large cosmological N-body simula- 
tion MICE (|Crocce et al.ll2010l ). They examined DM halos 
at 2=0.0 fc 0.5, searching for a halo pair matching the ini- 
tial conditions of lMastropietro fc Burkerd (|2008l ). They con- 
cluded that ACDM is excluded by more than 99.91% confi- 
dence level at z—0. Their results at 2=0.5 are inconclusive 
due to limited statistics. However, by fitting their pairwise 
velocity probability distribution function to a Gaussian dis- 
tribution, they were able to estimate the probability of find- 
ing a pair with W12 > 3000 km s"'^ to be 3.6 x 10~^ and 
?;i2 > 2000 km s~^ to be 2.2 x 10"^ at 2=0.5. They did 
warn that one must be careful about this approach since 
they are probing the tail of the distribution where their fits 
may not be accurate. 

Most recently iForero-Romero et al.l (|201[| ) approached 
the problem from a diff'erent pers pective. They studied data 
from the MareNostrum Universe (|Gottlober fc Yepesll2007l ) 
which contains baryonic matter in addition to collisionless 
DM. Instead of examining the pairwise velocities of DM halo 
pairs, they concerned themselves with the physical separa- 
tion between the dominant gas clump and its predominant 
DM structure. They argued that their approach provides a 
more robust comparison to observation; deriving the rela- 
tive velocity from the observations includes statistical and 
systematic uncertainties whereas the separation uncertainty 
is dominated by statistical errors in the measuring process. 
Additionally they point out that current simulations do not 
include the proper prescriptions for cooling, star formation, 
or feedback which implies that their predictions of the de- 
tailed X-ray properties of hot gas in massive halos are not 
robust. Using their method they found that large displace- 
ments between gas & DM are common in ACDM simulations 
therefore, 1E0657-56 should not be considered a challenge. 

In this pap e r, we take a similar approach to that of 
iLee fc Komatsul (|2010f ). and examine large ACDM A''-body 
simulations to see how common these high relative veloci- 
ties are among massive DM halos. One of the things that the 
earlier works have not performed is an examination of reso- 
lution and box size effect. Therefore we first conduct a study 
to determine the effects of increasing resolution or varying 
box sizes on the parameters of interest. We then examine 
our largest simulation in se arch for a pair matching t he ini - 
tial conditions required by JMastropietro fc Burkerd (120081 ) 
to reproduce the observed properties of 1E0657-56. 

The rest of the paper is organized as follows: Section[2] 
discusses simulation parameters, Section[3] shows the simu- 
lation results and examines the distribution of parameters 
relevant to this study, such as halo masses, pairwise velocity, 
and pair separation distances. Section|4] examines the simu- 
lation results at earlier redshifts of 2=0.296 fc 2=0.489, and 
how they relate to the bullet system. Finally, Section[5] con- 
tains concluding remarks and discussion of future prospects. 
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2 SIMULATIONS 

For our simulation s we use the GADGET-3 code (origi- 
nally described in ISpringell |2005|) which simulates large 
A''-body problems by means of calculating gravitational 
interactions with a hierarchi cal multipole expansion. It 
uses a particle-mesh method (IHocknev fc Eastwood! Il98ll : 
iKlvpin fc ShandarinllQS^ : IWhite et al.lll983h For ioiig-range 
forces and a tree method ( Barnes fc Hud 119861 ) for short- 
range forces. 

Cosmological parameters consistent with the cosmolog- 
ical constraints from the Wilkinson Microwave An isotropy 
Probe (WMAP) data and an iEisenstein fc Hul (|l998l ) transfer 
function were employed when creating the initial conditions 
for each simulation with random Gaus sian phases: (ilm, ^a , 
gp, -8, ns) = {Q.2&, 0.74, 72, 0.8, 1.0) (|Komatsu et al.ll2009l . 
[201l|). We note that we used a value of ns = 1.0 although the 
best-fit value from the WMAP data is ns=0.96. Two addi- 
tional simulations were ran with ns=0.96 and no differences 
in their high mass halo pairs were found from the ns = l sim- 
ulations, because the tilt in the primordial power spectrum 
mostly changes the small scale structures. All simulations 
contain only collision-less dark matter particles that inter- 
act solely through gravity. 

Several simulations with varying particle counts and 
box sizes were ran from z=100 to z=0. The list of simu- 
lations along with other parameters can be found in Ta- 
ble [T] Starting with the L250N125 run, the box size and 
particle count were simultaneously increased (from Lbox = 
250/i"^ Mpc to 2016/i"^ Mpc, and from iV = 125^ to 1008^ 
particles) in order to maintain the same mass resolution and 
gravitational softening length up until the L2016N1008 run. 
The second set of simulations were ran to examine the res- 
olution effect. We started with the original L250N125 sim- 
ulation and increased the particle count and decreased the 
gravitational softening length while keeping the box size the 
same, up to the L250N500 run. 



3 DATA ANALYSIS & RESULTS 

3.1 Halo Mass Function 

DM particles were grouped using a simplified version of 
the parallel friends-o f-friends (FOF) group finder SUBFIND 
ISpringel et al.M200lh . The code groups the particles into 
DM halos if they lie within a specified linking length (FOF 
LL). This linking length is a fraction of the initial mean 
inter-particle separation, for which we adopt a standard 
value of b=0.2. In order to be considered a halo it must 
also contain at least 32 particles. 

Figures [T] and [2] show DM halo mass functions in 
our si mulations. Both figures include the ISheth fc TormenI 
1 19991 ) mass f unction (ST) plotted as a black dotted line. Re- 
cent work by I More et al.l (|201in found that the commonly 
used value of 6=0.2 selects a significantly larger local over- 
density (5fof) than previously thought. Normally it is as- 
sumed that fe=0.2 results in (5fof ~ 60 (corresponding to the 
enclosed overdensity of 5 ~ 180), but their study finds that 
it results in 5fof ~ 80.61 which is a ~35% increase. We find 
that our mass function is slightly higher than the ST mass 
function on all mass scales. By regrouping the L1000N500 
sim using 6=0.1 we under-predict the number density on all 
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Figure 1. DM halo mass function at ^=0. This figure shows 
the box size effect; how increasing the simulation box size allows 
for a larger number of high mass halos. The abscissa uses a bin 
size of AlogMhai o=0.1. The black dotte d line is the ST mass 
function using the lEisenstein &: Hj l ll99a) transfer function. The 
solid magenta line is from the lGpcN500 simulation grouped with 
a linking length parameter of 6=0.1 instead of 0.2. 
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Figure 2. DM halo mass function at z=0. This figure shows the 
resolution effect; how increasing the resolution of a simulation 
allows for a greater number of small mass halos. The abscissa uses 
a bin size of A log Mt,aio=0-1. The black do tted line is the ST mass 
function using the lEisenstein &: Hul l ll99a) transfer function. The 
solid magenta line is from the 250MpcN500 run, using a linking 
length of fe=0.1 instead of 0.2. 
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Figure 3. Pairwisc velocity function at z = 0, demonstrating 
tlie box size effect. Each panel contains three lines representing 
pair separation distances of di2 = 2, 5, & 10 Mpc. Increasing the 
box size allows a higher «;i2 for pairs within di2 < 10 Mpc, while 
pairs within di2 < 5 Mpc only see a minor increase. Pairs residing 
within di2 < 2 Mpc see the smallest increase in D12 as the box 
size increases. 



Figure 4. Pairwisc velocity function at z=0, demonstrating the 
resolution effect. Each panel contains three lines representing pair 
separation distances of di2 = 2, 5, & 10 Mpc. Each subsequent 
increase in resolution allows for smaller structures to be resolved, 
leading to the increase in i'i2 at all separation distances. 



mass scales, as shown by the solid magenta line in Figures [T] 
&[2l Changes in b certainly have a significant impact on the 
halo mass function. 

Figure [T] shows that the number of high mass halos 
increases by increasing the box size from 250 h~ Mpc to 
2016 h~^Mpc while maintaining the same resolution. The 
lowest mass halo in all simulations shown in Figure [T] is 
Afhaiomin=l-84 X 10^''^ h"'^ Mq . The run with the largest box 
size (L2016) shows a slight shortage in the number of low 
mass halos around Mhaio=i lO^^'^* - 10"-^''/i"^Mo when 
compared to the other three runs with smaller box sizes. 
The most likely explanation for this shortage is that the 
lower mass halos were absorbed into higher mass halos. 

Higher resolution runs can resolve larger number of 
low mass halos as seen in Figure (2] The least massive 
halo for the highest resolution simulation (L250N500) has 
.^^haiomin=2.87 X IO^^H'^Mq, wMch is roughly two orders of 
magnitude lower than the lowest mass halos found in Fig- 
ure [T] 

While searching for a bullet-like pair of halos with 



masses on the order of Mbuiiet 



M, 



parent ; 



Figures [T] & [2] 



indicate that it is possible to form such massive halos in box 
sizes as small as 250 h~^Mpc at z=0 but there will be a low 
number of them. 



3.2 Pairwise Velocity Function 

In this section, we present the results on the pairwise velocity 
(^12 ~ \vi — V2\) function, i.e., the number of halo pairs 
within a velocity bin per unit volume {dn/dvi2). Figures [3] 
& [4] show dn/dvi2 with four panels for different simulation 



runs, each panel containing three lines for halo pairs with a 
separation distance of less than di2 = 2, 5, & 10 Mpc. 

Figure [3] shows that increasing the box size with a fixed 
resolution allows for a greater number of high i;i2 pairs, 
but with greater separation distances. Doubling the box 
size from L250 to L500 yields only a small increase in high 
Wi2 pairs. Doubling it again to LIOOO gives us a consider- 
able jump in high t;i2 pairs with separation distances of 
5< di2 <10 Mpc, while the 2< di2 <5 Mpc range only 
sees a moderate increase. Doubling the box size one final 
time to L2016, we again only see a moderate increase in 1112 
similar to going from the L250 to L500 sim. The number 
of close halo pairs with di2 <2 Mpc seem to remain fairly 
constant with relatively low W12 throughout changes in the 
box size. This implies that increasing the box size does not 
increase ?;i2 for pairs within 2 Mpc of one another. 

By increasing the resolution, the number of halo pairs 
with high U12 increases (Figure |3| , but unlike the case of 
enlarging the box, this does not necessarily come at the cost 
of increased separation distances. Each increase in resolu- 
tion gives us a larger number of low and high «i2 pairs on 
all distance scales. When compared to Figure |3] the sim- 
ulations shown in Figure |4] are better at resolving smaller 
structures and length scales, leading to larger values of U12. 
Unfortunately this data does not give us any information 
on the mass of the halos pairs, so increasing the resolution 
in order to increase the number of close high V12 pairs may 
not be beneficial when searching for high mass pairs such as 
1E0657-56. 
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Figure 5. Pairwise velocity vs. average mass of DM halo pairs at 
2=0. Here we show the box size effect; increasing the simulation 
box size increases the number of low-mass, high-j;i2 pairs more 
than the high-mass, high-tii2 pairs. Each increase in the box size 
and particle count yields better statistics, broadening the distri- 
bution of V\2- 



Figure 6. Pairwise velocity vs. average DM halo pair mass at 
2=0. This illustrates the resolution effect; how increasing the res- 
olution probes lower mass halo pairs. There is a slight increase in 
high-mass, high-j;i2 pairs, but the majority of the increase is in 
the low mass halos. As the particle count increases we can resolve 
smaller structures with higher Vi2- 



3.3 Relative Halo Velocity & Halo Mass 

It is useful to study the effects of different box sizes and res- 
olutions on the average mass of a halo pair vs. v\2- Figure [5] 
shows how increasing the box size with a constant resolution 
increases the number of low-mass, high v\2 halo pairs, along 
with increasing the number of high-mass, high-ui2 pairs to 
a lesser degree. As the box size increases, we are allowing 
for a greater number of rare high v\2 halo pairs which probe 
the tail of the distribution. 

Figure [6] shows that an increase in the resolution results 
in a larger number of low-mass, high-ui2 pairs, and a less 
substantial increase in the number of high-mass, high-i;i2 
pairs. Increasing the box size yields high vi2 pairs with in- 
creasing mass, while increasing the resolution yields a larger 
number of high v\2 pairs at the maximum halo mass allowed 
by the box. 



3.4 Cumulative i;i2 Function 

To examine how the box size and resolution affects the ac- 
tual number of high v\2 halo pairs, we plot the cumulative 
vi2 distribution function as shown in Figures [7] & (8] Chang- 
ing the box size (Figure (Tj) extends the curve to higher v\2- 
The larger box and particle count result in better statistics, 
which allows us to probe the high velocity tail of the v\2 
distribution as mentioned in the previous section. 

By increasing the resolution alone (Figure [8]), we see 
that the normalization of the cumulative «i2 distribution 
function becomes higher due to larger number of lower mass 
halos. These figures suggest that by increasing the box size 
and/or resolution one would be able to produce a halo pair 
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Figure 7. Cumulative Vi2 function of DM halos at 2=0. This 
figure shows how increasing the box size increases the number of 
high-t)i2 pairs, extending the tail of the distribution. 
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Figure 8. Cumulative j;i2 function of DM halos at z=0. This 
figure siiows tfie resolution effect. As tlie resolution increases, tiie 
normalization of the distribution increases due to a larger number 
of lower mass halos with higher velocities. 



Figure 9. The average halo peculiar velocity of five simulations 
used in this study, compared with the normalized prediction of 
linear theory described by Eq. JTJ. The data agrees with theory 
at 2 > 1, but the velocities begin to level off at z < 1, which is 
likely due to the virialization of the halos. 



with a greater wi2, however, as we saw earlier in Figures [5] & 
[S] the majority of high-i'i2 pairs have lower average masses 
than 1E0657-56. 



4 RESULTS AT EARLIER REDSHIFTS 

To be fully consistent with the observations of 1E0657-56, 
comparing our simulations at the same redshift as 1E0657- 
56 would be ideal. Up until this point, we have examined 
only simulation data at z—0, yet 1E0657-56 is observed at 
2=0.296. This difference in time of ~3.31 billion years can 
have a considerable impact on the velocities, sizes, and sepa- 
ration distances of the DM halos contained in the simulation. 
Another problem arises when we consider how we group the 
DM particles. At 2=0.296 the separation between the two 
halos of 1E0657-56 is di2 — 0.72 Mpc, which is larger than 
the linking lengths listed in Table[T](0.1-0.4 Mpc) for each of 
our simulations. At first glance it may appear that we could 
identify each halo independently within our sims, but when 
one considers their large masses, we find that this is not 
the case. The virial radius of each halo is found to be 1.42 
& 3.06 Mpc for the 'bullet' (Mbuiiet i; 1.5 x IO^Mq) and 
its 'parent' (Mparent — 1.5 x IO^^Mq), respectively. When 
two halos of this size are separated by ~0.72 Mpc, they will 
easily overlap, resulting in the EOF group finder identifying 
them as a single halo at the observed redshift of 2=0.296. 
If we assume the separation distan ce of 5 Mpc and infall ve- 
locity of 3000 km s~^ as required bv lMastropietro fc BurkertI 
(|200g ) to reproduce the observed quantities of 1E0657-56, 
then a halo pair in this initial configuration should be found 
at 2=0.489. 



4.1 Peculiar Velocities 

Before we examine the simulation at 2=0.489, we first com- 
pare our simulations to the prediction of linear theory for 
further validation. Linear theory predicts that for an Ein- 
stein de-Sitter (EdS) universe the growing mode of the pe- 
culiar velocity field grows as t^'^. The pecul iar velocity o f 
each mode in a non-EdS universe is given by (|Peebleslll98Gr ) 



H{z) a^ dP 
47r da 



(1) 



where H{z)=HoE{z) is the Hubble parameter, a is the scale 
factor, D is the growth factor for linear perturbations, and 

£(2) = [Q„,0(1 + zf + {1- Qkfi - f^m,o)(l + zf + nA,o]l/^ 

The peculiar velocity of each halo in five of our runs was 
calculated and averaged up to 2=10, then compared against 
the normalized theory curve in Figure [5] Our simulations 
agree well between 2=6 to 2=1.0, but start to deviate from 
the linear theory curve at 2 < 1.0, which is likely due to 
their virialization. 



4.2 Pairwise Velocity: Linear Theory 

Ijuszkiewicz et al.l (|l999r ) proposed a simple closed-form ex- 
pression relating the mean relative velocity of pairs of galax- 
ies at a fixed separation to the two-point correlation function 
of mass density fluctuations: 



gf»|/l|l + o, 



(2) 



where H is the Hubble parameter, r = ax is the proper 
separation, / = dlnD/dlna, a « 1.2 — 0.657, 7 is the log- 
arithmic slope of the correlation function, ^ = C/[1 + C]j 
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Figure 10. Auto-corrclation function of DM iialos in the 
L250N500 run at 2=0—6. The vertical cyan dashed hnes indi- 
cate r=l, 3, 5, & 10 Mpc, where we measure the evolution of § 
as a function of redshift. Symbols lying along these dashed lines 
represent the ^-values used in Eq. ^ for producing the dashed 
lines in Figure [TTI For comparison, we also show the dashed black 
line with a slope of g ex r~^'^ — the result consistent with the 
z=0 SDSS galaxies jZehavi et al.ll2010l) . 



^ = 'ix~^ j^ ^y^dy, and ^ is the two point correlation func- 
tion. At 2=0, the value of / is ~0.5, and then it asymptotes 
to unity at 2:>8. 

To obtain theoretical results based on Eq. ^ that can 
be compared with our simulations, we calculate ^ by corre- 
lating the center-of- mass positions of halos with a random 
data set and use the lLandv fc Szalavl (|l993l ) estimator 



C(r).: 



DD - 2DR + RR 
RR ' 



(3) 



where DD, DR, & RR represents halo pair counts for Data- 
Data, Data-Random, & Random-Random data sets at a 
given value of r. The result of ^(r)haio for the L250N500 sim 
is plotted in Figure [TO] Higher values of ^haio correspond to a 
larger probability that another halo lies at a separation of r. 
The value of ^haio decreases with increasing r, implying that 
halos tend to cluster more on smaller scales. The value of 
Chaio also decreases with increasing redshift, meaning halos 
are less clustered in the earlier universe. 

To compare our simulation with Eq. ((2]), we calculated 
the average pairwise halo velocities (W12) for pairs residing 
within physical shells of IMpc thickness (±0.5 Mpc) around 
r=l, 3, 5, & 10 Mpc for the L250N500 run. The results are 
shown in Figure [TTI where the solid curves represent simu- 
lation data, the dashed curves correspond to the theoretical 
predictions of Eq. ((2]) using ^-values from Figure 1101 and 
th e different colors disting uish between different values of 
r. Ijuszkiewicz et al.l (|l999l ) did not consider the effect of 
galaxy bias relative to dark matter, and without any cor- 
rection, we find that (1)12) of halos in our simulation are 



Figure 11. Solid lines: Average pairwise halo velocities {v\2) 
from the 250MpcN500 run residing in physical shells of IMpc 
thickness with the indicated radii. Dashed lines: Theoretical {v-12) 
curves given by Equation (|2]l using the f values from Figure 1101 
at each co rresponding radius . The dashed cyan line represents 
data from iFukushige fc Sutol 11200 Ih at a separation distance of 
r=1.52Mpc. When these curves reside below unity the Hubble 
flow is greater than their pairwise velocities, at unity their phys- 
ical separations remain constant, and above unity their pairwise 
velocities arc greater than the Hubble flow. 



somewhat higher than those predicted by Eq. ([2]). Therefore 
we invoke an ad hoc correction factor of x 1.5 to account for 
this effect, and the dashed lines in Figure [TT] include this 
multiplication factor in the right-hand-size of Eq. ((2)1 . After 
this correction, our simulation agrees with Eq. ^ very well 
for r = 3 & 5 Mpc, but there is some deviation from theory 
for the r = 1 & 10 Mpc results. The shape of the theory curve 
is determined by the competition between increasing H{z), 
decreasing ^, and increasing / with increasing redshift. 

IFukushige fc Sutol (|200ll ) examined the validity and 
limitations of the stable condition {—V12/ Hr = 1), which 
states that the mean physical separation r of galaxy pairs 
is constant on small scales. They found a significant time 
variation in the mean pairwise peculiar velocities and ar- 
gued that this behavior was not due to a numerical artifact, 
but a natural consequence of the continuous merging pro- 
cess. This irregular oscillatory behavior could be reduced by 
averaging over cosmological volumes larger than 200 Mpc'', 
resulting in a more accurate estimate of the mean pairwise 
veloci ty. Our data is also consistent with IFukushige fc Sutol 
(|200H ) (dashed cyan line in Figure [TT) in that the oscillatory 
behavior is suppressed due to our cosmological volume being 
greater than 200 Mpc"', and their result for r=1.52Mpc lies 
between our r=l & 3 Mpc curves. 
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Figure 12. Pairwisc velocity function for the L2016N1008 run 
at ^=0.0, 0.296, & 0.489. There is a slight increase in the number 
of pairs at the highest end of the i;i2 distribution as the redshift 



increases. 



4.3 In Search of the 'Bullet' 

Hereafter we will only be examining our largest simulation 
(L2016N1008) at redshifts of 2=0.0, 0.296, and 0.489. In Fig- 
ure [I2l we show the redshift evolution of the pairwise veloc- 
ity function {dn/dv\2) from z=Q to 2=0.489. Qualitatively 
this plot changes very little with redshift, except that there 
is a slight increase in the number of pairs at the highest end 
of the Wi2 distribution. Pairs within separation distances of 
rfi2 < 2 Mpc have maximum v\2 on the order of ~1800 km 
s~^ at 2=0.296 and 2=0.489. For pairs with greater di2, the 
maximum Wi2 reaches as high as ~3300 km s^'^. 

In Figure 1131 we show the redshift evolution of the av- 
erage halo mass vs. their pairwise velocity. One can see the 
effect of halo mergers, and the number of high-mass halo 
pairs with (Mhaio) > 10^^ M© are increasing from 2=0.489 
to 2=0. The cyan dashed lines in the 2=0.489 panel illus- 
iss o 

., ity c 

by iMastropietro fc BurkertI (|2008l '). Two pairs are found in 
our simulation near the region of interest, but their masses 
and velocities are still too low. 



4-3.1 Candidate Halo Pairs 

Table [2] hsts the five halo pairs with highest (i\/haio) for 
2=0, 2=0.296, and 2=0.489. A simulation of this size 
(comoving 2/i^^Gpc) produces many halo pairs massive 
enough to match that of 1E0657-56 at the examined red- 
shifts. While the separation distances of these pairs may 
be in the range we are interested in, the pairwise veloci- 
ties are too low to match the required Di2=3000kms^^ by 
IMastropietro fc BurkertI (|2008h . 

Table [3] lists the five halo pairs with the highest Vi2 



Figure 13. Average mass of halo pairs vs. their pairwrise velocity 
for the L2016N1008 run at 2=0.0, 0.296, & 0.489. In the bottom 
panel (^=0.489) the horizontal dashed line represents an average 
pair mass of 8.25 X lO^^'M© for 1E0657-56, and the vertical dashed 
line represents a pairwi s e velo city of 3000 km s~^ suggested by 
IMastropietro fc Burkei^ bOOSi) . 



at the three examined redshifts. All halo pairs in this table 
match or exceed the required ui2 of 3000 km s~^, but they 
miss the mark when it comes to the other observables of 
1E0657-56. All of the halos in this table have masses one 
or two orders of magnitude lower than Mbuiict & A/parent. 
The mass ratios are also a bit high; the lowest being ~0.3 at 
2=0.489 compared to 0.1 for 1E0657-56 at 2=0.296. None 
of the collision angles are head on, yet most are highly in- 
clined. Lastly the separation distance of cadi pair at 2=0.489 
is somewhat large; jMastropictro & Burkcrt (2008|) set their 
initial separation at proper 5 Mpc while each pair in this 
table is separated by >7.5 Mpc. 



4-3.2 Simulation Requirements to Produce the 'Bullet' 

In Figures [7] & [51 we examined the cumulative v\2 distribu- 
tion, however these figures included a large number of low- 
mass halos which are of little interest to this study. There- 
fore in Figure 1141 we restrict the halo sample to those with 
masses greater than W^Mq at 2=0, 2=0.296, & 2=0.489. 
With increasing redshift we see a decrease in the total num- 
ber density of halo pairs above IQ^'^Mq. 

Assuming that the trend of the cumulative v\2 func- 
tion would continue to higher velocities with increasing 
box size (as was the case for 2=0 shown in Figure [7]), we 
can fit a line to the 2=0.489 curve and estimate the box 
size and particle count required to produce at least one 
halo pair with a specified ui2. A quadratic of the form 
y = yo + as + bx^ was fit to the 2=0.489 curve between 
the values of i;i2=800— 1500kms~^, and we obtain the best 



fit values of 2/o=-3.97, a=-3.31 x 10" 



b=1.59 X 10-^ 
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Figure 14. Comoving number density of lialo pairs in the N2016N1008 run with masses above IO^^Mq at 2=0, 0.296,& 0.489. We 
also over-plot a quadratic fit described in the text for z=0.489. The horizontal dashed hne illustrates the number density of halos writh 
v\2 = 3000 km s~^ corresponding to a box size of (4.48/i~^ Gpc)'^ and 2240'^ DM particles. The black filled circles represent the 1112 
values listed in Table |4] 



Based on this fit, we estimate the minimuin box sizes and 
particle counts (for the same resolution as the L2016N1008 
run) required to prod uce at least one halo w i th th e ini- 
tial velocities given by Mastropietro fc Burkerd (|200a ) , and 
ISpringel fc FarraJ (|2007l '). The result is listed in Table g] 

Our result suggests that we would need a simulation 
box size of (4.48/i"^Gpc)^ & 2240^ DM particles in order to 
produce at least one halo pair with an average mass greater 
than IO^Mq and V12 > 3000 kms"^ at z=0.489. The exact 
values of the required box size and particle count is some- 
what sensitive to the range of U12 used for the fit, therefore 
the values listed in Table 3] should be taken as a rough esti- 
mate. The required simulations are so large and they would 
take significant computational resources which is currently 
not feasible for us. 



4-3.3 Probability of Finding the 'Bullet' 

We also examine the probability distribution function 
(PDF) of «i2 for halos with (A/haio) > lO^^M©. We per- 
form a least square fit to the data usin g a skewed normal 
distribution (jAzzalini fc Capitanidl2009l ). and calculate the 
probability of finding a halo pair with W12 > 3000 km s^^ at 
z=0.489. 



In Figure [TSl we show the binned PDF data with blue 
circles, and the best-fit skew normal distribution as the red 
curve. By integrating the PDF from D12 = 3000 km s~^ to 
infinity, we calculate the probability of finding a halo pair 
with masses greater than lO'^* Af© and 1112 > 3000 km s~^ 
to be P(> 3000 km s~^) = 2.8 x 10"*, which is roughly 
one order of magnitu de higher than calculations done by 
iLee fc Komatsul (|201(]| ') (P=3.6xl0"^). This very low prob- 
ability corroborates our earlier finding that it is very dif- 
ficult to produce a massive halo pair with a high 1)12 
matching the required i nitial configuration suggested by 
iMastropietro fc Burkerd (|2008l ). 



5 CONCLUSIONS 

We performed many A''-Body cosmological simulations with 
varying box sizes and resolutions in order to examine how 
changing these parameters affect the search for high-ui2 halo 
pairs comparable to the initial conditions required to repro- 
duce the observed properties of the 1E0657-56 system in 
non-cosmological simulations. Using our largest L2016N1008 
run, we examined the pairwise velocities, halo masses, and 
halo separation distances at 2=0.0, 0.296, fc 0.489. 

We find that the high- 1112 tail of the distribution extends 
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to a greater velocities as we increase the simulation box size. 
We also find that the number of high-t)i2 pairs increased as 
we increase the particle count and resolution with a fixed box 
size, however, this increase is mostly due to lower mass halos 
which do not correspond to the characteristics of 1E0657-56. 
We find that the redshift evolution effect is not very strong 
for the Vi2 distribution function. 

As we show in Table [31 some of the halo pairs have 
a high relative velocity similar to the initial conditions re- 
quired to reproduce the observational quantities of 1E0657- 
56 in non-cosmological simulations, but they are galaxy 
group-scale halos (10^^ — lO^^Mo) and much less massive 
than the observed estimates for 1E0657-56. 

We find that, in A^-body simulations with comoving 
volumes of less than {2h~^ Gpc)^, it is very difficult to re- 
produce a system that resembles the initial conditions re- 
quired to reproduce the observational properties of 1E0657- 
56. Based on the extrapolation of our cumulative v\2 func- 
tion, we find that one needs a simulation with a comoving 
box size of (4.48 h~^ Gpc)^ and 2240^ DM particles in order 
to produce at least one pai r of halos that resembles t he in i- 
tial conditions suggested bv lMastropietro fc Burkerd (|2008l ). 
In the future it would be useful to run larger simulations 
(e.g., with ~5Gpc box and ^2500'^ particles) to improve 
the statistics of massive halos. 

From the simulated V\2 PDF of halos, we calculated the 
probability of finding a halo pair with iii2 J? 3000 km s^^ 
and masses J5 IQ^'^Mq to be 2.76 x 10~^, which is som ewhat 
larger than previous work bv lLee fc Komatsul ((20101). How- 
ever, both probabilities are quite small and the difference is 
negligible. These results suggest that a system like 1E0657- 
56 is currently incompatible with the concordance ACDM 
universe, if its initial condition really re quires an initial 
pair w ise velocity of vi2 ^ 3000 km s~^. As lLee fc Komatsul 
(|2010l ) discussed in detail, there seems to be more systems 
like 1E0657-56 being observed already, which exacerbates 
the incompatibility in terms of probability. One other pos- 
sibility is that there is something wrong with the referred 
non-cosmological simulations, and the suggested initial t;i2 
must be revised to a lower value. 




Figure 15. Pairwisc velocity probability distribution function 
for halo pairs with masses above 10^'^Mq in our L2016N1008 
run. The blue circles represent i)i2 binned PDF data, the blue 
curve is the linearly interpolated va lues, and the red curve is th e 
best-fit skew normal distribution I IAzzalini fc Capitaniol 120091) . 
Integrating the fit from Vi2 = 3000 km s"'^ to infinity gives 
P(> 3000kms-i) = 2.8 X 10"*. This very low probability sug- 
gests that it is very difficult to produce a halo pair with high mass 
and high-j;i2 as the observed 1E0657-56. 
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Table 1. Summary of Simulations 



Run Name 



Box Size 

[ft-i Mpc] 



Particle Count 



-1 



[h~^ Mq] [h-^ kpc] 



FOF LL 
[h-^ kpc] 



lox Size Effects 




L250 N125 


250 


L500 N250 


500 


LIOOO N500 


1000 


L2016 N1008 


2016 


Resolution Effects 




L250 N125 


250 


L250 N165 


250 


L250 N250 


250 


L250 N500 


250 



125^ 
250^ 
500^ 
1008^ 



125^ 
165^ 
250^ 
500=* 



5.74 X 10" 


80 


400 


5.74 X 10" 


80 


400 


5.74 X 10" 


80 


400 


5.74 X 10" 


80 


400 


5.74 X 10" 


80 


400 


2.50 X 10" 


60.6 


303 


7.18 X 10"' 


40 


200 


8.97 X 10^ 


20 


100 



Note. — Summary of simulations employed in this paper. M^m is the mass of each DM 
particle, e is the comoving gravitational softening length, and FOF LL is the friends-of- 
friends linking length. The top four simulations explore the effects of increasing box size 
with fixed resolution, while the bottom four explore the effects of increasing resolution 
with a fixed box size. 



Table 2. Highest Mass Pairs 



Pair 


1'12 


e 


A/i 


A/2 


Mass Ratio 


d 


ri virial 


I'2 virial 


2=0 

1 


1670 


165 


5.71E+15 


5.02E+14 


0.088 


8.70 


5.67 


2.52 


2 


1792 


46 


5.71E+15 


1.99E+14 


0.035 


7.84 


5.67 


1.85 


3 


1767 


75 


5.71E+15 


l.OlE+14 


0.018 


7.63 


5.67 


1.48 


4 


1624 


80 


5.71E+15 


7.33E+13 


0.013 


7.13 


5.67 


1.33 


5 


2316 


72 


5.71E+15 


7.04E+13 


0.012 


6.20 


5.67 


1.31 


2=0.296 


















6 


1360 


141 


3.80E+15 


3.50E+14 


0.092 


9.55 


4.18 


1.89 


7 


1533 


44 


3.80E+15 


2.61E+14 


0.069 


6.23 


4.18 


1.71 


8 


1486 


56 


3.80E+15 


2.51E+14 


0.066 


10.00 


4.18 


1.69 


9 


1425 


129 


3.80E+15 


2.13E+14 


0.056 


6.20 


4.18 


1.60 


10 


2007 


112 


3.80E+15 


1.78E+14 


0.047 


5.65 


4.18 


1.51 


2=0.489 


















11 


869 


91 


3.28E+15 


5.59E+14 


0.170 


8.78 


3.70 


2.05 


12 


1277 


111 


2.64E+15 


1.07E+15 


0.405 


8.11 


3.44 


2.55 


13 


1875 


132 


2.45E+15 


1.19E+15 


0.485 


3.86 


3.36 


2.64 


14 


1257 


108 


2.45E+15 


1.08E+15 


0.440 


4.83 


3.36 


2.55 


15 


1256 


54 


3.28E+15 


1.73E+14 


0.053 


6.01 


3.70 


1.39 



Note. — Five halo pairs with the highest average halo mass from the L2016N1008 simula- 
tion at 2=0.0, 2=0.296 and 2=0.489. The values of v\2 are given in kms~^, collision angles 
9 in degrees, masses (Mi,M2) in Mq, pair separation distances (^12) and virial radius of 
each halo in physical Mpc. Although this simulation can produce massive pairs matching 
the observed mass of 1E0657-56, these pairs have too low relative velocities, and too large 
separation distances. 
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Table 3. Highest Velocity Pairs 



Pair 


Vi2 


e 


A/i 


M2 


Mass Ratio 


d 


^1 virial 


r2 virial 


z=0 
31 


3674 


103 


3.64E+13 


2.71E+13 


0.746 


8.83 


1.05 


0.95 


32 


3199 


151 


2.14E+13 


2.02E+13 


0.946 


8.20 


0.88 


0.86 


33 


3133 


134 


5.83E+13 


2.60E+13 


0.446 


9.09 


1.23 


0.94 


34 


3095 


113 


8.20E+13 


4.56E+13 


0.556 


9.21 


1.38 


1.13 


35 


3053 


108 


8.20E+13 


2.14E+13 


0.261 


9.11 


1.38 


0.88 


z=0.296 


















36 


3538 


143 


3.35E+13 


1.96E+13 


0.586 


9.94 


0.86 


0.72 


37 


3282 


125 


4.96E+13 


2.37E+13 


0.477 


9.39 


0.98 


0.77 


38 


3141 


155 


8.60E+13 


3.41E+13 


0.396 


8.80 


1.18 


0.87 


39 


3089 


170 


6.93E+13 


2.77E+13 


0.400 


5.27 


1.10 


0.81 


40 


3053 


153 


4.16E+13 


2.48E+13 


0.597 


8.60 


0.93 


0.78 


2=0.489 


















41 


3361 


128 


6.75E+13 


2.60E+13 


0.385 


8.81 


1.01 


0.74 


42 


3312 


148 


5.66E+13 


3.18E+13 


0.561 


8.03 


0.96 


0.79 


43 


3239 


102 


6.75E+13 


2.37E+13 


0.350 


7.57 


1.01 


0.72 


44 


3109 


146 


2.94E+13 


2.37E+13 


0.804 


9.57 


0.77 


0.72 


45 


3083 


103 


7.56E+13 


2.37E+13 


0.313 


9.25 


1.05 


0.72 



Note. — Five halo pairs with highest v\2 found in the L2016N1008 simulation at 2=0.0, 
2=0.296 and 2=0.489. The values of v\2 are given in km s~^, collision angles 8 in degrees, 
masses (Mi, M2) in Mq, pair separation distances ((^12) and viral radius of each halo in physical 
Mpc. None of these high velocity halo pairs are massive enough to match the observations of 
1E0657-56. 



Table 4. Simulation Requirements to produce a Bullet 



Reference 


^"12 
[km s-i] 


Box Size 

[/i-i Mpc] 


Particle Count 


Mastropietro & Burkert ("2008) 
Sorinffel & Farrar ('2007) 


3000 
2057 


4480 
2224 


2240^ 
1112^ 



Note. — Required box size and particle number needed to produce at least 
one halo pair with an average mass greater than IO^^Mq and a certain value 
of v\2 at 2=0.489 suggested by each of the authors. See text in ii l4.3.2l for more 
details. 
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